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Abstract 

We discuss a simulation algorithm for dynamical fermions, which combines the multi- 
boson technique with the Hybrid Monte Carlo algorithm. The algorithm turns out to give 
a substantial gain over standard methods in practical simulations and to be suitable for 
dealing with fermion zero modes in a clean and controllable way. 

1 Introduction 

In this contribution, we discuss a simulation algorithm, which we call Polynomial Hybrid Monte 
Carlo (PHMC) for full lattice QCD, implemented in the case of Wilson fermions. The al- 
gorithm is based on two key ingredients, the interplay of which appears to be crucial for a per- 
formance gain over the standard HMC algorithm. The generation of gauge-field configurations 
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is performed as first proposed in 0]: the usual inverse of the Dirac operator is approximated 
(as in Liischer's multiboson algorithm 0) by a polynomial in the operator and, with this poly- 
nomial defining the fermion action, the simulation is done using a standard small-step Monte 
Carlo method. The second crucial ingredient is the way of correcting for the above-mentioned 
polynomial approximation. We suggest to do this by means of an efficient reweighing tech- 
nique, which is reminescent of earlier ideas, such as [Q, and allows us to deal in a clean and 
controllable way with fermion zero modes, as discussed below. 



2 The PHMC algorithm 



Denoting by U the configuration of lattice gauge links, the expectation value of any gauge 
invariant observable O = 0[U], in full QCD with nj = 2 degenerate fiavours, may be written 

as 

(1) 



(C) = Z-^ J VUe~^^^^^det{Q^[U])0[U] 



where Sg is the standard plaquette action for the pure gauge sector (/3 = G/g^) and Q is the 
Dirac operator for Wilson fermions multiplied by 75: 



Q[U]x,y = Co75 



1 + -^f^cswO'^uFfj_J^ 5x,y 



(2) 



with = 1 ±7/^, Ffj,u the standard 'clover' discretization of F^^,, and = Q. Here we present 
simulation results only for csw = 0. For csw > the computational cost per single molecular 
dynamics trajectory has been evaluated and performance tests are in progress. 

We use the polynomial approximation of described in the polynomial of degree 

n, denoted as Pn^e{s), approximates in the range e < s < 1 with a relative fit error bounded 
by 5 = 2 (where e > 0). Choosing cq such that \\Q^\\ < 1, we may write the 

corresponding polynomial of in a factorized form: 

n 

PnAQ')=Cn,el[iQ-^l)iQ'^'^) ^ (3) 

k=l 

where Cn,e is a positive constant, = y/zk = fik + it^k and the ZkS are the complex roots of 
Pn^s), as in ||^. For the values of n and e used in our tests, a careful ordering of the monomial 
factors appearing in (^ was essential, at least on computers with 32-bit precision, in order to 
keep to a negligible level the rounding error in constructing Pn^Q"^) using the factorized form, 
eq.(i. 

The full QCD [uf = 2) partition function may be represented as 

Z = j VUV(p^V(pVri^VT]We-^^^+^P^ 

Sp = Sp[u,<p,v] = <P^PnAQ'[u])<P + v^v (4) 



by introducing the auxiliary pseudofermion fields (i.e. boson fields with spin and colour indices) 

0, rj and the 'correction' factor W = Wit], U]: 

W = exp {r^\l - IQ' . P^,iQ')]-')v} . (5) 

In the PHMC algorithm, the update of the gauge field configuration is performed by using 
the ^-version ^ of the HMC algorithm for the 'approximate' but still non-local action Sg + Sp 
of eq. (I). We denote averages in the theory with action Sg + Sp as (. . .)p: reweighing with W 
yields the true averages, denoted as (. . .), for any observable O = 0[U]: 

{0) = {W)-p'{OW)p. (6) 

The number (Ncorr) of evaluations of W per single molecular dynamics trajectory (updating U) 
is relevant for the level of statistical error on (O), although (O) itself is correct, within statistical 
uncertainties, for any value of N^orr- Each evaluation of W requires a trivial Gaussian update 
of the ?7-field and the solution of the system [Q^P„^e(Q^)]x = V- 

3 The results 

All the results discussed here, for both HMC and PHMC algorithms, have been obtained using 
the even-odd preconditioned form of the Q operator, denoted by and a Sexton- Weingarten 
leap-frog integration scheme. We adopted Schrodinger functional boundary conditions and 
monitored few observables: the plaquette (P), the lowest (Amm) and the highest {\max) eigen- 
values of (normalized in such a way that Xmax ~ !)• 

Several tests on the 4^ lattice with different values of (ra, e) showed that the chosen polyno- 
mial approximation should not be too bad, in order to avoid that reweighing observables with 
W induce large statistical fluctuations in true averages. The approximation should not be too 
good either, in order to keep as low as possible the computational cost per single trajectory 
(of length ~ 1). In practice, a reasonable compromise is obtained by choosing e ~ 2(Amm) 
and n such that 5 ~ 0.01-0.02, which means n scaling as e~^/^. This criterion (plus short runs 
monitoring the statistical fluctuations of W) allows us to quickly choose reasonable values of 
(n,e). 

In Table 1, we compare HMC and PHMC algorithms, as far as the computational cost 
per single trajectory and the relative statistical error for P and Xmin are concerned. Bare 
lattice parameters are (3 = 6.4, n = 0.15 on the 4^ lattice. On the 8^ lattice we choose 
(3 = 5.6 and k = 0.1585 ^ Kc, which are the same bare parameters used for comparing the 
multiboson technique and the HMC algorithm in ref. @]. The chosen values of {n,e, Ncorr) 
are (12, 0.036, 1) for the 4'^ lattice, corresponding to one of the best choices, and (48, 0.026, 2) 
for the 8^ lattice. The values of e reflect a factor larger than 10 for the condition number of 
Q^, between the 8^ and 4^ lattices. Further details concerning the results reported in Table 

1, such as the proper definition of the computational costs, molecular dynamics parameters 
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Algorithm 


Cq0 


a{P)/{P) 




44 


PHMC 


540 


0.00024 


0.0064 


4^ 


HMC 


868 


0.00020 


0.0057 


8^ 


PHMC 


3974 


0.00027 


0.037 


8^ 


HMC 


7398 


0.00021 


0.039 



Table 1: We give the single trajectory cost, Cq^, in units of Q(f) operations. We compare the 
relative statistical errors for P and Xmin, obtained from a jack-knife blocking analysis, with 
18000 and 2745 trajectories for 4'^ and 8'^ lattices, respectively. Note that we estimate a 10 to 
20% uncertainty on the statistical errors. 



and average values with statistical uncertainties, may be found in ref. |I|]. The formulae for 
computational costs appearing there are also valid for the 0(a)-improved fermions {csw > 0). 

Looking mainly at the 8'^ lattice case, which is the one with a relatively large condition 
number of about 700, our results for the PHMC algorithm may be summarized as follows. 
First, average values of measured observables always agree with the corresponding ones from 
the HMC algorithm (within statistical uncertainties); without reweighing with W, systematic 
deviations, increasing with 1 — {W)p, are observed, as expected. Second, for suitable choices 
of (n, e), the statistical errors of measured observables are the same as for the HMC algorithm, 
at least within their relative uncertainty, which we estimate to be about 10-20%. At the same 
time, the computational cost per single trajectory, Cq^, is almost a factor of 2 lower than for 
the HMC algorithm. This is a consequence of the fact that Cq^ is basically proportional to 
for the PHMC and to A~|„ for the HMC algorithm. Last but not least, we find that PHMC 
and HMC algorithms sample configuration space in a somewhat different way: for instance, the 
low-lying end of the distributions of Xmin look different (Fig. 1). 

As expected from the properties of the polynomial approximation of (Q^)"^, the PHMC 
algorithm generates configurations with low-lying modes {Xmin < of (Q^) with a higher prob- 
ability than algorithms like HMC or the multiboson technique made exact by an accept/reject 
step 0: in particular the PHMC algorithm generates, with small, non-zero probability, con- 
figurations carrying fermion zero modes, which occur with vanishing probability (i.e. never in 
a finite run time) when using the other algorithms mentioned above. However, fermion zero 
modes give a finite, non-zero contribution to all those observables where the divergence of quark 
propagators compensates for the vanishing of the probability measure. In ref. 0], we discuss 
how we may in principle deal with fermion zero modes, in order to evaluate the correction 
factor W and quark propagators, by making use of existing minimization algorithms (tested 
up to Xmin ~ Xmax) ■ Therefore, the PHMC algorithm looks particularly suitable to study 
the contribution of low-lying fermion modes to physical observables in full QCD. For the same 
reasons, one may also expect that the PHMC algorithm may overcome energy barriers, related 
to low-lying fermion modes, easier than the HMC algorithm, leading to a better exploration of 
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Figure 1: The distribution of A = XminiQ'^) as obtained with equal statistics on the 8^ lattice 
for both PHMC and HMC algorithms. 

configuration space. 

This work is part of the ALPHA collaboration research programme. We thank DESY for 
allocating computer time to this project. 



References 



[1] 

[2] 

[3] 
[4] 
[5] 
[6] 
[7] 



R. Frezzotti and K. Jansen, Phys. Lett. B402 (1997) 328. 

Ph. de Forcrand and T. Takaishi, Nucl. Phys. B (Proc. Suppl.) 53 (1997) 968. 

M. Liischer, Nucl. Phys. B418 (1994) 637. 

B. Bunk et. al., Nucl. Phys. B (Proc. Suppl.) 42 (1995) 49. 

S. Gottlieb et. al., Phys. Rev. D35 (1987) 2531. 

K. Jansen, Nucl. Phys. B (Proc. Suppl.) 53 (1997) 127. 

A. Borrelh, Ph. de Forcrand and A. Galh, Nucl. Phys. B477 (1996) 809. 



5 



